This notebook calculates the acoustic complexity index (ACI) as defined in Pieretti, et al. 2011. The ACI is based on the "observation that many biotic sounds, such as bird songs, are char- acterized by an intrinsic variability of intensities, while some types of human generated noise (such as car passing or airplane transit) present very constant intensity values. (Pieretti, et al. 2011)."
Pieretti, N., A. Farina, and D. Morri. 2011. A new methodology to infer the singing activity of an avian community: The Acoustic Complexity Index (ACI). Ecological Indicators 11: 868-873. doi: 10.1016/j.ecolind.2010.11.005
In [1]:
import numpy as np
from numba import guvectorize, float64
from nacoustik import Wave
from nacoustik.noise import remove_background_noise
from nacoustik.spectrum import psd
import matplotlib.pyplot as plt
%matplotlib inline
In [2]:
filepath = "/Users/Jake/Desktop/seewave_examples/160317-145000_24k.wav"
In [3]:
sound = Wave(filepath)
sound.read()
sound.normalize()
# trim sound to first 20 seconds
sound.samples = sound.samples[0:480000, :]
sound.n_samples = 480000
f, t, a = psd(sound)
In [4]:
# figure configuration
dpi = 192
channels = sound.n_channels
fig, ax = plt.subplots(channels, 1)
fig.set_dpi(dpi)
fig.set_figwidth((920 / dpi) * 3)
fig.set_figheight((460 / dpi) * 3)
plt.subplots_adjust(left=0, bottom=0, right=1, top=1, wspace=0, hspace=0)
fig.set_frameon(False)
# specify frequency bins (width of 1 kiloherz)
bins = np.arange(0, (sound.rate / 2), 1000)
# plot spectrogram for each channel
for channel in range(channels):
spec = ax[channel].pcolormesh(t, f, a[channel], cmap='viridis')
ax[channel].set(ylim=([0, sound.rate / 2]),
yticks = bins.astype(np.int) + 1000)
ax[channel].tick_params(length=12, color='white',
bottom=False, labelbottom=False,
top=False, labeltop=False,
labelleft=False,
labelright=False)
ax[channel].set_frame_on(False)
returns the ACI for each frequency band for each channel
implemented as a universal function via numba.guvectorize
In [5]:
@guvectorize([(float64[:,:,:], float64[:], float64[:], float64[:,:])],
'(c,f,t),(),()->(c,f)', nopython=True)
def calculate_aci(a, time_delta, block_duration, aci):
block_delta = int(np.around(block_duration[0] / time_delta[0]))
n_blocks = int(np.floor(a.shape[2] / block_delta))
remainder = int(a.shape[2] - (block_delta * n_blocks))
for channel in range(a.shape[0]):
for f_band in range(a.shape[1]):
aci_f_band = np.empty(shape=(n_blocks + 1))
for block in range(n_blocks):
d = np.empty(shape=(block_delta - 1))
for t_step in range(block_delta - 1):
d[t_step] = np.abs(a[channel, f_band, (t_step * (block + 1))] - \
a[channel, f_band, ((t_step * (block + 1)) + 1)])
D = d.sum()
aci_f_band[block] = D / a[channel, f_band, \
(block_delta * block):(block_delta * (block + 1))].sum()
if remainder > 1:
d = np.empty(shape=(remainder - 1))
for t_step in range(remainder - 1):
d[t_step] = np.abs(a[channel, f_band, ((t_step * n_blocks) + t_step)] - \
a[channel, f_band, ((t_step * n_blocks) + t_step + 1)])
D = d.sum()
aci_f_band[-1] = D / a[channel, f_band, -(remainder + 1):-1].sum()
# average aci value over blocks
aci[channel, f_band] = aci_f_band.sum() / (n_blocks + (remainder / block_delta))
convert the spectrogram values to watts before computing the ACI
In [6]:
aci = calculate_aci(10**(a / 10),
(t[1] - t[0]),
5.,
np.empty(shape=(a.shape[0], a.shape[1])))
In [7]:
colors = ['red', 'blue']
fig, ax = plt.subplots(figsize=(15, 5))
for channel in range(aci.shape[0]):
p = plt.plot(f, aci[channel], color=colors[channel],
label="channel {0}".format(channel))
title = ax.set_title('ACI')
xlabel = ax.set_xlabel('frequency (herz)')
ylabel = ax.set_ylabel('ACI value (dimensionless)')
legend = ax.legend()
ax.grid(True, which='major')
remove background noise
In [8]:
ale = remove_background_noise(a, iterations=5)
plot spectrogram after noise removal
In [9]:
# figure configuration
dpi = 192
channels = sound.n_channels
fig, ax = plt.subplots(channels, 1)
fig.set_dpi(dpi)
fig.set_figwidth((920 / dpi) * 3)
fig.set_figheight((460 / dpi) * 3)
plt.subplots_adjust(left=0, bottom=0, right=1, top=1, wspace=0, hspace=0)
fig.set_frameon(False)
# specify frequency bins (width of 1 kiloherz)
bins = np.arange(0, (sound.rate / 2), 1000)
# plot spectrogram for each channel
for channel in range(channels):
spec = ax[channel].pcolormesh(t, f, ale[channel], cmap='viridis')
ax[channel].set(ylim=([0, sound.rate / 2]),
yticks = bins.astype(np.int) + 1000)
ax[channel].tick_params(length=12, color='white',
bottom=False, labelbottom=False,
top=False, labeltop=False,
labelleft=False,
labelright=False)
ax[channel].set_frame_on(False)
compute ACI
In [10]:
aci_ale = calculate_aci(10**(ale / 10),
(t[1] - t[0]),
5.,
np.empty(shape=(ale.shape[0], ale.shape[1])))
plot the ACI
In [11]:
colors = ['red', 'blue']
fig, ax = plt.subplots(figsize=(15, 5))
for channel in range(aci.shape[0]):
p = plt.plot(f, aci_ale[channel], color=colors[channel],
label="channel {0}".format(channel))
title = ax.set_title('ACI (after noise removal)')
xlabel = ax.set_xlabel('frequency (herz)')
ylabel = ax.set_ylabel('ACI value (dimensionless)')
legend = ax.legend()
ax.grid(True, which='major')